; docformat = 'rst'
;
;+
;
; :Purpose:
;   Output lifetime information of optimized model
;   Estimates steady state stratospheric, tropospheric and OH lifetimes
;
; :Outputs:
;   Output stored in <A2I_PATH>/output/lifetimes.csv
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, May 7, 2013
;
;-
pro a2i_lifetime, parameters, lifetime, model

  yearToSec=365.25*24.*3600.

  
  stratospheric_lifetime=fltarr(parameters['NPOLLUTANTS'])
  tropospheric_lifetime=fltarr(parameters['NPOLLUTANTS'])
  OH_lifetime=fltarr(parameters['NPOLLUTANTS'])
  global_lifetime=fltarr(parameters['NPOLLUTANTS'])
  
  
  for pi=0, parameters['NPOLLUTANTS']-1 do begin

    ;nYEARS of 10, or less (if sim is less than 10 years)
    ; HOWEVER, we remove the last year below so it's actually nyears-1
    nYears=min([10, n_elements(model['LIFETIME', pi, 0, *])/4])
    nIt=10
    mLifetime=fltarr(12, nIt*(nYEars-1)*4)
    mOH=fltarr(8, nIt*(nYEars-1)*12)
    mEmissions=fltarr(4, nIt*(nYEars-1)*12)
    mTscale=fltarr(17, nIt*(nYEars-1)*4)
    
    for it=0, nIt-1 do begin
      mLifetime[*, it*(nYears-1)*4:(it+1)*(nYears-1)*4-1]=model['LIFETIME', pi, *, -4*nYears:-5]
      mOH[*, it*(nYears-1)*12:(it+1)*(nYears-1)*12-1]=model['OH_SCALE', *, -12*nYears:-13]
      mEmissions[*, it*(nYears-1)*12:(it+1)*(nYears-1)*12-1]=model['EMISSIONS', pi, *, -12*nYears:-13]
    endfor

    if parameters['TRANSPORT_ESTIMATE'] then begin
      print, 'WARNING! A2I_LIFETIME: add transport parameters'
    endif

    ;Run model
    conc=agage_box_model(model['IC', pi], mEmissions, $
      arr_OH=[parameters['OH_A', pi],parameters['OH_ER', pi]], mol_mass=parameters['MOL_MASS', pi], $
      lifeTime=mLifetime, $
      oh_scale=mOH, $
      overall_lifetime=lifetime_global, $
      transport_input_Extension=parameters['TRANSPORT_FILE'], $
      OH_input_extension=parameters['OH_FILE'])
    conc=agage_box_model(model['IC', pi], mEmissions, $
      arr_OH=[parameters['OH_A', pi],parameters['OH_ER', pi]], mol_mass=parameters['MOL_MASS', pi], $
      lifeTime=mLifetime, $
      oh_scale=mOH, $
      overall_lifetime=lifetime_strat, output_lifetime_boxes=[8, 9, 10, 11], $
      transport_input_Extension=parameters['TRANSPORT_FILE'], $
      OH_input_extension=parameters['OH_FILE'])
    conc=agage_box_model(model['IC', pi], mEmissions, $
      arr_OH=[parameters['OH_A', pi],parameters['OH_ER', pi]], mol_mass=parameters['MOL_MASS', pi], $
      lifeTime=mLifetime, $
      oh_scale=mOH, $
      overall_lifetime=lifetime_trop, output_lifetime_boxes=indgen(8), $
      transport_input_Extension=parameters['TRANSPORT_FILE'], $
      OH_input_extension=parameters['OH_FILE'])
    conc=agage_box_model(model['IC', pi], mEmissions, $
      arr_OH=[parameters['OH_A', pi],parameters['OH_ER', pi]], mol_mass=parameters['MOL_MASS', pi], $
      lifeTime=mLifetime, $
      oh_scale=mOH, $
      overall_lifetime=lifetime_OH, /output_oh_lifetime, $
      transport_input_Extension=parameters['TRANSPORT_FILE'], $
      OH_input_extension=parameters['OH_FILE'])

    global_lifetime[pi]=mean(1./(1./lifetime_global[-120:-1] + 1./lifetime_OH[-120:-1]))/365.25/24./3600.
    stratospheric_lifetime[pi]=mean(lifetime_strat[-120:-1])/365.25/24./3600.
    OH_lifetime[pi]=mean(lifetime_oh[-120:-1])/365.25/24./3600.
    tropospheric_lifetime[pi]=mean(lifetime_trop[-120:-1])/365.25/24./3600.

  endfor
  
  
  wh=where(finite(global_lifetime) eq 0 or global_lifetime gt 1.e6 or global_lifetime le 0., count)
  if count gt 0 then begin
    global_lifetime[wh]=1.e6
  endif
  wh=where(finite(OH_lifetime) eq 0 or OH_lifetime gt 1.e6 or OH_lifetime le 0., count)
  if count gt 0 then begin
    OH_lifetime[wh]=1.e6
  endif
  wh=where(finite(tropospheric_lifetime) eq 0 or tropospheric_lifetime gt 1.e6 or tropospheric_lifetime le 0., count)
  if count gt 0 then begin
    tropospheric_lifetime[wh]=1.e6
  endif
  wh=where(finite(stratospheric_lifetime) eq 0 or stratospheric_lifetime gt 1.e6 or stratospheric_lifetime le 0., count)
  if count gt 0 then begin
    stratospheric_lifetime[wh]=1.e6
  endif
  
  openw, fun, a2i_filestr(/Input, parameters['CASE_NAME'] +  '/output/lifetimes.csv'), /get_lun
  
  printf, fun, "Species, Transient, Global, OH, Stratospheric, Tropospheric,"
  
  for pi=0, parameters['NPOLLUTANTS']-1 do begin
    wh=where(lifetime['POLLUTANT'] eq pi and lifetime['MODEL'] gt 0., count)
    if count gt 0 then begin
;      meanLife=mean(1./lifetime['MODEL', wh])/yearToSec
      printf, fun, parameters['POLLUTANTS', pi], 0., global_lifetime[pi], oh_lifetime[pi], $
        stratospheric_lifetime[pi], tropospheric_lifetime[pi], format="(A, ',', 5(E16.8, ','))"
      print, strCompress(parameters['POLLUTANTS', pi] + ' lifetime: ' + $
        string(global_lifetime[pi], format='(f7.1)' ) $
        + ' years')
    endif else begin
      printf, fun, parameters['POLLUTANTS', pi], ",,,,,,"
      print, strCompress(parameters['POLLUTANTS', pi] + ' lifetime: ')
    endelse
  endfor

  close, fun
  free_lun, fun

end